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In this paper we discuss stability properties of various discretizations for axisymmetric systems 
including the so called cartoon method which was proposed by Alcubierre, Brandt et.al. for the 
simulation of such systems on Cartesian grids. We show that within the context of the method of 
lines such discretizations tend to be unstable unless one takes care in the way individual singular 
terms are treated. Examples are given for the linear axisymmetric wave equation in flat space. 
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I. INTRODUCTION 

Axisymmetric systems are notorious for the difficul- 
ties they pose in numerical simulations when they are 
expressed in coordinates adapted to the symmetry. The 
problem is due to the singular nature of the coordinates 
on the axis, i.e., the set of fixed points for the symmetry 
transformations or, equivalently, the set where the in- 
finitesimal generator of the symmetry - the Killing vec- 
tor - vanishes. Adapted coordinates, i.e., an angle (j) 
along the Killing vector, a radius r measuring 'distance' 
from the axis and a third coordinate z, yield coordinates 
(r, z) on the space of orbits of the action of the symme- 
try group. The general theory of group actions ^ 
yields the following facts. The orbit space contains two 
types of orbits jljl: the 'regular orbit' is a circle around 
the axis. It has maximal dimension (one) and the union 
of all regular orbits is dense in the space of orbits. The 
other type of orbit is a fixed point which has dimension 
zero. The set of fixed points is a submanifold of codimen- 
sion two, the 'axis'. The drop in dimension of the orbits 
shows up in the topology of the space of orbits which ac- 
quires a boundary corresponding to the axis points and 
in a degeneracy of the adapted coordinates: since on the 
axis (for fixed z) all values of address the same point 
the angular coordinate becomes irrelevant on the axis, so 
that the coordinate system degenerates there. 

Invariant tensor fields transform in a rigid way (Lie 
transport) under the symmetry so that they are deter- 
mined on an entire regular orbit once they are known 
on one of its points. Therefore, in order to determine 
the fields on the entire space it is enough to determine 
them on a hypersurface which is transversal to the or- 
bits. Thus, the problem is reduced to a problem in a 
space with lower dimension. For invariant tensor fields 
at a fixed point the invariance implies a transformation 
law among the tensor components which must be inter- 
preted as a regularity condition for the tensor field on the 
axis. 

The advantage of using the adapted coordinates is, of 
course, the manifestation of this reduction. In these coor- 
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dinates the components of the tensor fields do not depend 
on the angular coordinate along the KiUing vector. Thus, 
there are only two essential coordinates left over and the 
problem reduces in complexity because an axisymmetric 
problem on some open 3-dimensional domain reduces to 
an equivalent 2D-problem on a domain with boundary. 
The boundary conditions are obtained from regularity 
conditions on the axis. 

The disadvantage of using the adapted coordinates is 
that the reduced equations become formally singular on 
the axis. The regularity conditions on the axis guarantee 
that the individual terms in the equation have in fact a 
regular Hmit on the axis so there is not a real problem. 
Yet, numerically, one has to evaluate terms which are 
formally 0/0 or c» — oo which does pose problems. Thus, 
it seems that one has the choice between a 3D system 
with regular, such as Cartesian, coordinates and regular 
equations or a 2D system with degenerate coordinates 
and singular equations. 

In an attempt to combine the best of both worlds, Al- 
cubierre et al. proposed what they called the cartoon 
method. This scheme was designed to borrow from the 
singularity-free nature of full 3D Cartesian coordinates 
X, y, z and to allow the treatment of axisymmetric prob- 
lems without the memory constraints of a full 3D prob- 
lem. The method uses the Killing transport equation 
along the orbits to find the field values on the grid points 
of a Cartesian grid from its values on the hypersurface 
y = 0. These values are obtained by interpolation from 
the grid points lying in that hypersurface. For a second 
order method one needs to keep only three hypersurfaces 
y = const in memory so that the problem scales quadrat- 
ically with the number of points and not cubically. 

The plan of the paper is as follows. In section II we 
discuss the cartoon method in more detail. In order to 
have a specific problem at hand we apply it to the ax- 
isymmetric linear wave equation in section III. We write 
this equation as a first order symmetric hyperboHc sys- 
tem. The reason for this apparent complication is that 
we have implemented this method for the axisymmetric 
conformal field equations which are a symmetric hyper- 
bolic system and because the Hnear wave equation served 
as a model equation to test various implementations [||. 
It turns out that the crucial ingredient to the method 
is the interpolation and we discuss various possibilities. 
Furthermore, we show that there is no real difference be- 
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tween the cartoon method using Cartesian coordinates 
and the formulation in adapted coordinates so that the 
issues which are relevant for the former apply to the latter 
as well. In section |^ we analyse the stability properties 
of the cartoon method with different interpolation pro- 
cedures as well as other possible discretizations when the 
method of Hues is used for time evolution. We end with 
a brief discussion of the results. 



respect to the action of 5*0(2). The infinitesimal version 
of this condition is the equation of KiUing transport 

L^T = 0. (3) 

Suppose we want to solve a quasi-linear system of PDEs 
of the form 

t = F{daT,T,t,^) (4) 



II. THE CARTOON METHOD 

In this section we present the cartoon method in more 
detail. We consider a 3-dimensional Euclidean space £^ 
with Cartesian coordinates x = (x, y, z) on which we de- 
fine the action of the circle group SO{2) = G 

G={9={l-:)\a^ + b' = l] 
in the usual way by 
$ : (G,f^) 

(g, x) 1-^ x) = [ax — by, bx + ay, z) = g ■ x. 

(1) 

The infinitesimal generator of this action is the vector 
field ^ = xdy — ydx which vanishes on the axis a; = y = 0, 
the set of fixed points for the action. The action $ is the 
prototype for any axisymmetric system in the neighbour- 
hood of the axis in the sense that for any of these systems 
one can find local coordinates in which the action takes 
the above form ||. Therefore, we will focus here on 
fiat space because we are mainly interested in the local 
properties in a neighbourhood of the axis. Curvature will 
not play an essential role. 

Introducing adapted cyHndrical coordinates (r, cj), z) 
the coordinate expression for the action is 

and the Killing vector is given by 

These coordinates are singular on the axis which can be 
seen e.g., from the fact that the Jacobi matrix of the 
transformation between cylindrical and Cartesian coor- 
dinates has vanishing determinant there. 

A tensor field T in £^ is axisymmetric if for any g G G 
it satisfies the condition 

$*r = T (2) 

where $* denotes the pull-back with : £^ ^ £^ ,x i-^ 
g ■ x. This condition expresses the invariance of T with 



where T is some given time dependent axisymmetric ten- 
sor field on Here, da denotes covariant differentiation. 
The usual way to treat this system is to express it in 
terms of the adapted coordinates and to use the invari- 
ance condition (||) which asserts, that the components 
of T in the basis of adapted coordinates do not depend 
on (j). Thus, the essential coordinates are r and z and 
the problem is reduced to a 2-dimensional one. However, 
from the form of the covariant derivative operator, i.e., 
the Christoffel symbols of the metric expressed in cyHn- 
drical coordinates it is clear, that the equation will, in 
general, contain singular terms. 

A different procedure is to employ any coordinate sys- 
tem which is regular in a neighbourhood of the axis 
and then to make explicit use of the invariance condi- 
tion (||) or its infinitesimal version (||) So we may 
choose Cartesian coordinates with the axis at x = y = 0. 
Let T be an arbitrary axisymmetric tensor field on 
From the invariance condition (||) we can find its val- 
ues at any point in £^ from its values on the half-plane 
'H = {y = 0,x > 0}. To illustrate this let us first assume 
that T is a scalar field. From the invariance condition we 
have 

r(x) = T{g ■ x) 

for any rotation matrix g and each point x = {x,y,z). 
Thus, in particular, choosing 

( X V \ 

-y X 1 

for X 7^ 0, y 7^ we have <&g(x) = (r, 0, z) with r = 
\/ + y-^ so that 

T(x) =r(r,0,z). (5) 

Similarly, if T = Ti dx + T2dy + T3 dz is a 1-form, we 
have 

$;T(x) = Ti(.9 • x) d{<f;x) + T2{g ■ x) d{<f;y) + T^{g- x) d($;z). 
Inserting the explicit form of the matrix g yields 
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$*r(x) ~Ti{g ■ x) {adx — bdy) + T2{g ■ x) (bdx + ady) + T3{g ■ x) dz 

= {aTi{g ■ x) + bT2{g ■ x)) dx + {aT2{g ■ x) - bTi{g ■ x)) dy + T^ig ■ x) dz 

I 



so that, in this case, the invariance condition implies the 
three equations 

Ti(x) = aTi{g ■ x) + bT2{g ■ x), 
T2(x) =-6Ti(.g-x)+ar2(3-x), 
r3(x)=T3(g-x). 



Choosing, as before, a = x/^/x^ + y^ and b = 
— y/ + y^, we get the relations 

Ti(x) = i (a;Ti(r,0,z)-yr2(r,0,z)), 
r 

T2(x)==i (yri(r,0,z) + a-T2(r,0,z)), (6) 
T3(x)-T3(r,0,z). 

Similar relations hold for tensor fields with different type. 
This shows that we can get information about the values 
at all points outside the axis from values at the points on 
H except for the axis. On the axis these relations imply 
regularity conditions for the tensor fields. E.g., in the 
case of a 1-form above we get ri(0, 0,z) = T2(0,0, z) = 
which follows by taking limits from various directions 
towards the axis points. 

How are we to use this analytical background in numer- 
ical applications? Following [|| we cover the half-plane 
7i by a regular Cartesian grid 

Xi = iAx, Zk ~ kAz, with z, fc S Z, i > 

which we think of as being a part of a full 3D cartesian 
grid with points given by Xi,yj,Zk with j/o = 0. Each 
component of a tensor field T yields a grid function Tijk- 
The axisymmetry impHes that we should be able to de- 
termine the functions Tij^ from their values on Ti, alone. 

In order to treat equations like above one needs 
to compute approximations to the spatial derivatives of 
the components. This is straightforward in the case of 
X- and z-derivatives. However, for the y-derivatives we 
need to find the values of the grid functions at points 
outside Ti.. The situation is illustrated in Fig. 0. It shows 
the view from the positive z-axis down onto the (x, y)- 
plane. The black dots indicate grid points in Ti, the point 
on the left with the cross being the axis. The open circles 
indicate points outside H. with y = ±Ay. The values of 
the fields at these points have to be computed from the 
invariance equation (||). Following the orbits through 
these outside points one can connect them with points 
in Ti, indicated by the open squares, and one can find 
the values at the outside points from the values at these 




FIG. 1: Illustration of the method (see text) 



points on H. Thus, e.g., applying (||) in this situation we 
obtain for a scalar 

T(x„ ±Ay, zu) = T{^x^ + Ay^, 0, z,-), 

while the corresponding transformation for a 1-form fol- 
lows from (^). Note, that the squares are on Ti but they 
are not, in general, grid points. Thus, having reduced 
the problem to H the remaining question is how to de- 
termine the values at these points. The natural choice is 
to use some form of interpolation as it is also suggested 
in However, it is not clear what kind of interpolation 
should be used and it is the purpose of this work to point 
out that not all interpolation schemes perform well. 

III. THE AXISYMMETRIC WAVE EQUATION 

To have a concrete example at hand we will focus 
now on the axisymmetric scalar wave equation in fiat 
space. It is clear from the considerations above that the 
z-dependence is irrelevant for this analysis. Therefore, 
we simplify even further by assuming that the field is in 
fact axisymmetric and z-independent. Then the equation 
is 

where (j) ~ (j){t, x, y) is axisymmetric. Introducing the 
1-form d(j) ~ udt + V dx + w dy we can derive from this 
equation a symmetric hyperbolic first order system 

ii^Vx+Wy, (7) 
V ^ Ux, (8) 

W = Uy (9) 

with the constraint Vy — Wx = 0. The function 4> can be 
recovered from a solution of this system by either evolv- 
ing = M or by solving = Vx + Wy at each time 
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step. Since 4>{t,x,y) = (j){t, \J 
we have, in particular, a;, 0) 



= and w(i,0,0) = 0. 



Furthermore, taking derivatives we obtain u-y(f,x,0) = 
and ^^(tjXjO) = 0. Hence, we can forget about (|9|) and 
the constraint, because these equations are identically 
satisfied on "H. So we are left with the system 



U = Vx + Wy 

i) ^ Urr.. 



(10) 
(11) 



Before we go into the more numerical details it is worth- 
while to write down the system in adapted (i.e., cylindri- 
cal) coordinates (r, 4>, z) 



Vr 



(12) 



which follows by expressing the Laplace operator in these 
coordinates, using the independence of z and and then 
introducing the first derivatives as new dependent vari- 
ables. 

We solve these equations numerically by applying the 
method of lines, i.e., we discretize the spatial derivatives 
to obtain a system of ODEs. This can then be solved with 
standard ODE solvers. We put the equations on a grid 
as indicated above. Now Ti. is in effect a 1-dimensional 
grid with points Xi = iAx and we have grid functions 
Ui{t) = M(i, zAa;,0), Vi{t) = w(t, zAa;,0). According to 
the idea of the cartoon method we compute the term Wy, 
which of course does not vanish on by following the 
orbits. From the transformation law (^ for a 1-form we 
get 



w{t,x,y) 



y 



v{t, ^x^+y^O). (13) 



The discretization of the term Wy using centered differ- 
ences to get a second order accurate approximation yields 



Wy{t,Xi,0) 



w{t, Xj, Ay) - w(t, Xj, -Ay) _ 

2Ay x^ 



where = \/ xf + Ay^ and = v{t,x„,0). This, then, 
yields the final system obtained by application of the car- 
toon method 



U = Vx + —, V ^ Ux. 



(14) 



It is instructive to compare this system with the sys- 
tem dl^ ) obtained in adapted coordinates. Expanding in 
powers of Ay we get 

w(x*) Vi 1 f dv . . Vi\ Ay 

X 



Hence, up to the accuracy Ay^ for which the derivation 
of ( |l^ ) is valid, the two systems agree. 

One can, in fact, say even more. The infinitesimal 
invariance condition (||) or, equivalently, taking the y- 
derivative of (O) and evaluating at y = yields 



Wy {t,X,0) 



v(t, X, 0) 



for all a; > 0. Hence, whenever we use the cartoon 
method to approximate Wy to a certain order in the grid 
spacing Ay we approximate v/xio the same order. Thus, 
there is no difference in that order of accuracy between 
the systems ( |l^ and (|l^). And, therefore, the discus- 
sions in the following section apply to both. 



IV. NUMERICAL ISSUES 

In this section we discuss the stability properties of 
various implementations of the systems (|l|) and (^. 
In the spirit of what was said in the previous section 
we regard the term ^ as an approximation to ^ which 
is not necessarily located on a grid point. Then, the 
main question is how to approximate this term on TL. We 
implement several possibilities and look at their stability 
properties. 

In order to have a well defined problem we have to 
worry about boundary conditions at the origin x — 
and at an outer boundary which we put aX x = 1. The 
conditions at the origin follow from the transformation 
properties of the fields under the symmetry. They are 

v{t, -x, 0) = -u(i, X, 0) ^ v{t, 0, 0) = 0, (15) 

du 

u{t, -X, 0) = u{t, X, 0) —(t, 0, 0) = 0. (16) 

ox 

The fact, that a; = is the fixed point of the symmetry is 
reflected in the simultaneous vanishing of v and x. This 
also implies that at the origin we have ^ = Vx so that at 
X = the time evolution of uq — u{t, 0, 0) is given by the 
equation uq = 2vx{t,0,0). 

Since we are interested in the properties of the scheme 
in the interior we implement some sort of periodic bound- 
ary conditions. Since we want to avoid a second origin 
at a: = 1 we imagine that our wave equation lives on a 
2-sphere and is symmetric under rotations of the sphere 
around its poles. Now we impose the condition that it 
be also symmetric under reflection across the equator. 
This discrete symmetry allows us to consider the equa- 
tion only along a meridian from the north pole down to 
the equator. The symmetry implies conditions for the 
fields at the equator which we put at x = 1: 

v{t, 1 - a;, 0) = ~v{t, 1 -f- a;, 0) v{t, 1, 0) = 0, (17) 

Qu 

u(t,l~x,0) = u(t,l + x,0) ^ — (t, 1,0) = 0. (18) 

ox 

In the context of the method of lines we write down a 
system of ODEs for the grid functions and Ui. This 
will be a linear system of the form / = A/ and we can 
find its stability properties by looking at the spectrum of 
the matrix A which we determine numerically. 

With these preparations out of the way we can now 
write down the first of our discretizations for the sys- 
tem dl2). This is the first one which comes to mind. 
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namely 



UN 



2vi 

IT' 

Vi+l 



2h 

VN-l 



vo = 0, 



2h 



(19) 



VN = 0. 



Here we have chosen h = 1/N where + 1 is the number 
of grid points so that Xi = ih. The matrix A which 
corresponds to this discretization has dimensions 2(7V + 
1) X 2(A^ + 1) and the form 





^ 


C 




vD 






where C is the matrix 



/ 



2 
1 

_i 

2 



1 

-1 



1 
2 

0/ 



and D is 



D 



/ 



V 



\ 



1 i 

2 " 2 

■ • 



0/ 



This system of ODEs is solved by a standard ODE solver, 
such as a Runge-Kutta scheme or a multistep scheme 
like the Adams-Bashforth schemes. We also consider the 
'iterated Crank-Nicholson' (ICN) scheme which has be- 
come popular among numerical relativists (see ||^ for an 
analysis of this scheme). Relevant for deciding about 
the stability of a difference approximation for ODEs is 
the region of absolute stability (see e.g. ||^, |^). This is 
a closed set fl in the complex plane which consists of 
those z = XAt £ C for which the method appHed to the 
equation ii = Xu with step size Ai produces bounded 
approximations. The stability regions for the schemes 
mentioned above are shown in appendix |^ and the sta- 
bility region(s) for ICN is determined in appendix |^. 

Thus, in the context of the method of lines, the sta- 
bility conditions can be obtained by determining the 
eigenvalues of the discretization matrix A and checking 
whether it is possible to scale it so that it fits entirely into 
the stability region of the ODE solver. The scale factor 
determines the maximal admissible time step. Strictly 
speaking, this consideration applies only to the Hnear 
case. For non-Hnear equations this check should be un- 
derstood as a 'rule of thumb' |B|. 
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FIG. 2: Spectrum of the 'straightforward' discretisation, 
scheme 2a 



The spectrum aN{hA) of hA is shown for N = 100 
in Fig. g. We point out some of its properties which are 
valid for all other discretizations below as well. The spec- 
trum is symmetric with respect to reflections about the 
real axis as well as about the imaginary axis. While the 
former symmetry is due to the reality of the equation, 
the latter is a consequence of the fact that the underly- 
ing equation is a second order wave equation which has 
in- and outgoing modes. Furthermore, the spectrum is 
concentrated mostly on the imaginary axis with a few 
'outliers'. It can be seen that the eigenvalues with the 
largest, in absolute value, real part remain unchanged 
when we change N which implies that the spectral ra- 
dius of A is proportional to N. 

An interesting feature of the spectrum are the two 
'handles' which are located near ±i. These, and the out- 
liers are due to the v/x term in the equation and, depend- 
ing on the particular discretization scheme, they may or 
may not be present as we will see below. Since the eigen- 
values of the corresponding continuous system are given 
by iiAfc , where Xk is the k-th zero of the Bessel function 
ji, all the eigenvalues with nonvanishing real part are 
spurious, i.e., they do not have an analogue in the con- 
tinuous problem. Thus, the corresponding eigenvectors 
are parasitic modes because they do not approximate a 
regular mode. 

When combined with a time integrator this spectrum 
and the stability region of the integrator decide about the 
stability of the overall scheme. The stability regions of 
some ODE solvers are sketched in appendices |^and |^. 
Some of these include an interval on the imaginary axis 
while others include only the origin. All the stability re- 
gions extend further into the left half plane than into the 
right. This implies that if a spectrum contains points 
in the right plane (i.e., modes with positive real part) 
then these modes tend to either unnecessarily diminish 
the time step because they have to be scaled into the 
stability region (if the stability region of ODE solver ex- 
tends into the right half plane int he first place) or else 
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they create instabilities. 

Another possibiHty to overcome instabihties is the ad- 
dition of numerical dissipation. This has the effect that 
the spectrum of A is shifted towards the left. However, 
the presence of the handles and of the outliers implies 
that the amount of dissipation necessary for stabilising 
the scheme is rather high compared to cases where the 
spectrum is located entirely on the imaginary axis. This 
means that the evolution cannot be followed accurately 
because ultimately the waves will be damped out. Thus, 
the appearance of handles and outliers seems to leave 
us with the choice between the Scylla of instability and 
Charybdis of inaccuracy. Therefore, we take the appear- 
ance of these modes as a sign for an inadequate scheme. 

The discretizations which we have looked at can be 
divided into three groups. Discretizations of the sys- 
tem ( |l^ obtained from the cartoon method of computing 
the transversal derivatives, 2nd order discretizations of 
the system ( |l2|) and three somewhat non-standard dis- 
cretizations of that system. In the list given below we 
only show the formula for the inner points of the grid. 
The values at the boundary points are obtained by us- 
ing the appropriate symmetry conditions. Furthermore, 
we omit the equation for Vi in those cases where its dis- 
cretization is the same as for the derivative term in the 
Ui equation. 

1. Using the cartoon approximation we discretize the 
system (|lj) to second order in h. The derivatives 
are taken to be centered in all the three different 
schemes, while we use different interpolations for 
the term 
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FIG. 3: Two (dots) and three (circles) point interpolation of 
the cartoon method, schemes la, lb 



(a) Two point interpolation (see Fig. y): 



2h 



hx^ 



(b) Three point interpolation (see Fig. ^): 



Vj+l - Vj-l 

2h 



where Pi-i{x) = (x - Xi){x - Xt+i)/2h'^, 
Pi{x) = -{x - Xi+i){x - Xt-i)/h'^ and 

Pt+l{x) = (.T - Xi){x - Xi_l)/2/7,2. 

(c) Centered two point interpolation (see Fig. ||): 

iH+i - v,^i , {xi+i - x^)vi-i + {x^ - Xi^i)vij^i 



2h 



It is obvious that the interpolation in schemes la 
and 16 is not useful because of the appearance of 
the handles and the outliers. This suggests that 
for higher order Lagrange interpolation one has to 
expect the same kind of behaviour. In fact, in our 
studies related to we find that e.g. 5th order La- 
grange interpolation is as unstable as interpolation 



FIG. 4: Centered two point interpolation for the cartoon 
method, scheme Ic 



using Chebyshev polynomials. However, the fact 
that scheme Ic yields a purely imaginary spectrum 
(up to round-off error) suggests that the use of the 
centered interpolation is advantageous. In view of 
the fact that the cartoon method provides an ap- 
proximation of the 1/x term in (|lj) we now look 
at this system directly. 

2. Three second order discretizations of the sys- 
tem (H). 

(a) The scheme discussed already above (see 
Fig. I) 



2h 



ih 



(b) A centered approximation for the numerator 
of the second term (see Fig. ||) 



2h 



2ih 
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FIG. 5: Scheme 2b 



(c) A staggered grid discretization based on the 
variables Ui = u{t, ih), fi+1/2 = v(t, ih + h/2) 
(see Fig. I) 

= h + 2^h ' = 




FIG. 7: A fourth order centered scheme, 3a 



(b) The next scheme makes use of further infor- 
mation coming from the theory of group ac- 
tions. While the radial coordinate r is not 
a smooth coordinate on the orbit space its 
square is smooth. Hence, it seems natural to 
express the variables as functions of the co- 
ordinate X = r'^ . Then the minor redefinition 
u{x) = 2u(r) and v{x) — rv(r) of the variables 
yields the (regular) system 

dtu = dxv, dtv = xdxu, 

which is discretized in the standard way (see 
Fig. I) 



2h 



-{Ut+l ~ Ui-l). 



FIG. 6: Scheme 2c 



Again, we find that discretizing in a centered and 
compatible way yields a purely imaginary spectrum 
(again up to round-off error). The staggering in 
scheme 2c is suggested by the structure of the equa- 
tions which is ultimately related to the meaning of 
the variables. While u is a time derivative located 
at the grid points w is a spatial derivative which 
should be located between the grid points. 



2) are 



3. The final three discretization schemes of 
somewhat different from the ones above: 



(a) This scheme takes the idea of centering to 
fourth order and again yields a purely imagi- 
nary spectrum (see Fig. ^ 
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12h 



+ 



FIG. 

f2 



8: Standard discretization for redefined system 



6ih 
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(c) The same system discretized using a staggered 
grid (see Fig. ||) 



H+l 



/2 = {i + l/2)(iti - Ui-i). 



o 

D OO, 



a 2D axisymmetric code to a full 3D code because the 
underlying grid structure and the equations need not be 
changed. 

Clearly, this investigation cannot deliver the ultimate 
solution to the problem. It should be regarded as pro- 
viding a guide to ehminate 'bad' discretizations and sug- 
gesting possible 'good' ones. For instance, we have im- 
plemented the schemes Ic and 3a in our axisymmetric 
code 1^ for solving the conformal field equations. While 
these schemes perform well for the axisymmetric wave 
equation they still tend to produce slowly growing insta- 
bilities in the case of the quasilinear conformal field equa- 
tions. This behaviour might be related to the fact that in 
most cases the real parts of the eigenvalues vanish only 
up to round-off error and that these small real parts start 
to grow due to non-linear interactions. It would be use- 
ful to have a scheme in which the spectrum is exactly on 
the imaginary axis and whose only non-vanishing eigen- 
values are the physical ones. However, such a scheme is 
probably impossible to find. 



FIG. 9: Staggered discretization for redefined system 
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V. CONCLUSION 

We discussed in this work several schemes for treating 
axisymmetric systems which are difficult to handle nu- 
merically due to the degenarcy of adapted coordinates 
and hence the formal singularity of the ensuing equa- 
tions. We applied various schemes to the axisymmetric 
wave equation written as a symmetric hyperbolic sys- 
tem for the derivatives and determined the spectrum of 
the ensuing matrices. As a consequence of the shape 
of the stability regions of common ODE solvers we take 
the appearance of handles and outHers described in sec- 
tion IV as an indication for a 'bad' discretization scheme. 



These eigenvalues correspond to grid modes which, when 
evolved with a standard ODE solver, will lead to insta- 
bilities. In order to avoid them one needs to be careful 
in discretizing the equations because the derivative term 
and the 1/x term must be treated in a compatible way. 
This is due to their common geometric origin. In the 
general case, where tensors of higher rank are involved 
different terms may have to be combined. Exactly which 
terms these are is indicated by the behaviour of the ten- 
sor components under the symmetry transformation. 

One particular method we discussed is the cartoon 
method. It turns out that this method is not really dif- 
ferent from a direct implementation of the axisymmetry 
in terms of the adapted coordinates. Here the problem is 
to make the interpolation compatible with the discretiza- 
tion in the y = half plane so that parasitic modes are 
excluded. However, the advantage in using the cartoon 
method is that it admits a quick and easy extension of 



APPENDIX A: STABILITY REGIONS 

We illustrate the procedure for determining the stabil- 
ity region of a particular time evolution scheme for an 
explicit Runge-Kutta scheme. Let u = f{t,u) be the 
(system of) ODE(s) which we intend to solve. An ex- 
plicit Runge-Kutta scheme with s 'stages' to compute an 
approximation ui of u{to + h) is formulated as follows: 



ki = f{to,uo), 

k2 = f{to + C2h,uo 

h = f{to + csh^uo 



kg 



/(io + Cs/l, Mo 

Mo + h{biki + 



ha2iki), 

hasiki + ha32k2), 

husiki + has2k2 H h /la^s-ifcs-i), 

■ • + bgk„). 



The particular scheme is characterized entirely by the co- 
efficients aik, bi and Ck which are conveniently presented 
as a tableau 





C2 
C3 


021 
031 


032 






Cs 


a si 


Os2 


Oss-1 






bi 


62 


bs-1 


bs 



The method is said to be of p-th order if for sufficiently 
smooth / the condition 

\\ui - u{to + h)\\ < ChP+^ 
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FIG. 10: Stability regions for explicit Runge-Kutta schemes 
with p — s 



holds. 

To get to the stabihty region one apphes the scheme 
to the specific model equation u = \u where A S C is 
arbitrary. Going through the procedure above it is clear 
that the final result is of the form 

ui ~ P{Xh)uQ 

where P{z) is a polynomial of degree s, the so called sta- 
bility function of the method. Since ||ui|| = \P{Xh)\ \\uo\\ 
we find that the method will produce bounded approxi- 
mations if and only if |P(A/i)| < 1. This motivates the 
definition of the stability region as the set 

n = {zeC\ \P{z)\ < 1}. 

In case of an implicit scheme the stability function, 
obtained in the same way by applying the scheme to 
the model equation, is not a polynomial but a ratio- 
nal function. In a similar way, the stability regions of 
linear multistep schemes are obtained. Here, the sim- 
plest one is the leapfrog scheme which updates m„+i from 
the values on the previous time levels m„ and m„_i as 
Un+i = Un-1 + 2hf{tQ + nh,Un). It is easy to see that 
the stability region for the leapfrog scheme is the part of 
the imaginary axis between — i and i. 

The stability regions for explicit Runge-Kutta schemes 
of order p = s for s — 1, 2, 3, 4 are shown in Fig. |l^. The 
stability regions for the second and third order Adams- 
Bashforth methods are shown in Fig. |ll| as an example for 
stability regions of multistep methods. For the stability 
regions of other multistep schemes we refer to | pO| , pl| . 
Qualitatively, these regions extend into the left half plane 
of the complex plane including possibly some part of the 
imaginary axis. For higher order methods the stability 
region may even extend into the right half plane. 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 



FIG. 11: Stabihty regions for second (thin fine) and third 
(thick line) order Adams-Bashforth methods 



APPENDIX B: THE ICN TIME EVOLUTION 
SCHEME 

The ICN scheme was described in in application to 
a one-dimensional advection equation. Viewing the ICN 
scheme as a time stepper to solve systems of ODEs u — 
f{t, u) one can deduce from this description the following 
algorithm to compute Mi « M(to -I- h) 

h = f{to + \h, uq + i/ifci), 
h = f{to + \h, Uq + \hk2), 

Ui = Uq + hkg. 

Cleary, this is an explicit Runge-Kutta scheme with the 
tableau 
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1 
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2 








1 





1 
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2 






1 
2 








1 

2 















1. 



Note, that k iterations correspond to s = fc + 1 stages 
of the corresponding Runge-Kutta scheme. Application 
of this scheme to the model equation u ^ Xu yields the 
stability function 

P.,.,..+=g(l)'.i+,44g(i)', 

fe=0 fe=0 

The fact that the first three terms agree with the Taylor 
expansion of exp(z), the stability function of the exact 
time evolution, shows that this scheme is second order 
accurate 0. The stability regions are obtained as the 
level set {Ps{z)Ps{z) = 1}. The regions for some values 
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FIG. 12: The stability regions for ICN(s) for s = 1, 2, s = 3, 4, s = 5, 6 and s = 7, i 



of s are shown in Fig. We see that the imaginary 
axis intersects the stabihty regions for s = 3,4, 7, 8, . . . 
in a finite interval while for the other values they have 
only the origin in common with the imaginary axis. To 
explain this somewhat unexpected result we look at the 
behaviour of the modulus of the stability function on the 
imaginary axis in a neighbourhood of the origin. Writing 



and inserting z = \y into Q{z) = Ps{z)Ps{z) we find 



P,{z) = l + zJ2 



k=0 



1 + 2- 



1 - 



1 - 



-(i2//2) 



s+l( 



-1) 



s+l\ 



(ij//2r+2(i + (-in + 2(y/2) 



2s+2 



1 + (y/2)2 



r 



This shows that for all s > 1 



1 - 



,s+2 



, s+1 

1 ^ yS+1 



s even 
s odd. 



The behaviour of the stability regions is determined by 
the coefficient of the lowest order term. These lead to the 
sequence 1, 1/4, -1/4, -1/16, 1/16, ... for s = 1, 2, . . . . 
The negative coefficients imply a decrease of Q{iy) along 
the imaginary axis close to y = and hence they corre- 
spond to the cases where the stability region contains a 
part of the imaginary axis. The other cases, when the 
coefficient of the dominating term is positive, are those 
when the stability region just touches the imaginary axis 
at the origin. In these cases, systems with eigenvalues on 
or to the right of the imaginary axis such as in partic- 
ular the advection equation and other hyperbolic equa- 
tions cannot be stably evolved with ICN. Thus, we come 
to the same conclusion as in ||6|, namely that ICN(A:) 
for fc = 2, 3, 6, 7, . . . is unstable for hyperbolic equations. 
Furthermore, since the method is always second order ac- 
curate regardless of how many iterations are performed, 
there is no use in performing more than two. In fact, in 



most practical cases it is probably more efficient to use a 
higher order Runge-Kutta method instead. 

Another interesting aspect of the ICN scheme is its 
relationship with the Crank-Nicholson method. This is 
a second order accurate, implicit time stepping scheme 
whose stability function is given by 



P{z) 



l-Hz/2 
1 - z/2 



and whose stability region is the entire left plane. Tak- 
ing the limit of Ps{z) for s ^ oo we find that in fact 
Ps{z) — > P{z). However, absolute convergence holds 
only in the circle \z\ < 2. This has the following im- 
plication: iterating the ICN scheme until convergence 
to produce approximations to the exact Crank-Nicholson 
method works only within that circle. Thus, this scheme 
(which might be termed 'the full ICN scheme') also has 
a stability barrier and hence a limitation of the time step 
in contrast to the exact CN scheme which is uncondition- 
ally stable. Thus, we conclude that the ICN scheme is 
not an efficient alternative to higher order Runge-Kutta 
or implicit schemes like CN. 
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